function rnd = mnormal_rnd (r, c, z0)

   # usage:  mnormal_rnd (r, c[, z0])
   #
   # Return an r by size(c) matrix of normally distributed random vectors with
   # covariance c and, optional, truncated at z=z0.

   global PAR
   if ~isfield(PAR, "par_opt")
      PAR.par_opt = {"UniformOutput", false, "VerboseLevel", 0} ;
   endif
   
   if !(isscalar (r) && (r > 0) && (r == round (r)))
      error ("mnormal_rnd:  r must be a positive integer");
   endif
   if (size(c,1) != size(c,2))
      error ("mnormal_rnd:  c must be a quadratic matrix");
   endif

   [g p] = chol(c) ;
   if (p != 0)
      warning("xds:xds", "mnormal_rnd: matrix not pd, attempting bend\n") ;
      [g p] = chol(bend(c)) ;
      if (p != 0)
	 save /tmp/tt.ob c
	 error("mnormal_rnd: matrix still not pd, saving to /tmp/tt.ob\n") ;
      endif 
   endif


   if (nargin < 3)
      x = stdnormal_rnd (r, size(c,1));
   else
      zs0 = z0 / g ;
      zs0 = mat2cell(zs0, 1, ones(1, columns(zs0))) ;
      x = parcellfun_(@(z0) gstdnormal_rnd(r, z0), zs0, PAR.par_opt{:}) ;
      x = cell2mat(x) ;
   endif
   
   rnd = x * g ;

endfunction


function z = gstdnormal_rnd (r, zs0)

   ## usage:  z = gstdnormal_rnd (r, zs0)
   ##
   ## generate r N(0,1) numbers below zs0

   z = [] ;
   while (length(z) < r)
      w = randn(r, 1) ;
      I = w < zs0 ;
      z = [z ; w(I)] ;
   endwhile 
   z = z(1:r) ;

endfunction
